Statistical Mechanics of the Anomalous Behavior of Tetrahedral Liquids 
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Tetrahedral liquids such as water and silica-melt show unusual thermodynamic behavior such as 
a density maximum and an increase in specific-heat when cooled to low temperatures. There is a 
debate in the literature whether these phenomena stem from a phase transition into a low-density 
and high-density liquid phases, which occur in the supercooled regime. Here we consider a model 
of tetrahedral liquids for which we construct a volume-constrained statistical mechanical theory 
which quantifies the local structure of the liquid. We compare the theory to molecular dynamics 
simulations and show that the theory can rationalize the simulations semi-quantitatively. We show 
that the anomalous density and specific heat behavior arise naturally from this theory without 
exhibiting a liquid-liquid phase-transition. We explain that this theory may or may not have a 
phase transition, depending on the volume and temperature dependence of the energy and entropy 
which are sensitive to small changes in the parameters of the model. 
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"Tetrahedral liquids" are composed of molecules which 
tend to form strong bonds with four neighbors in a 
tetrahedral geometry. Several materials share this type 
of interaction, most notably water, silicon and silica. 
These liquids exhibit non-trivial thermodynamics at low 
temperature: when cooled, instead of steadily reducing 
in volume, they have a density maximum; the specific 
heat, usually a result of thermal fluctuations, increases 
rather than decreases and the thermal expansion becomes 
smaller and even reaches negative values. There is a 
debate in the literature whether these phenomena stem 
from a first-order phase transition into a low-density and 
high-density liquid phases or a continuous change fl. Ex- 
tensive studies used experiments [H Q , simulations [H, 0] 
and theory to clarify this point. Several groups have 
invoked different model descriptions which incorporate 
Hamiltonian terms that were believed to represent im- 
portant physical effects, cf. Refs 0, @ just to mention 
a few. A more comprehensive approach was used by 
Truskett at al. Q which used a statistical mechanical 
model to describe the properties of the two-dimensional 
Mercedez-Benz model. Despite of these efforts the origin 
of the anomalies is not completely clear. 

Here we apply statistical mechanics methods to study 
the issue on the basis of the three-dimentional Stillinger- 
Weber model Hamiltonian, which is a generic, empirical 
model of tetrahedral liquids. The Stillinger- Weber model 
incorporates two and three-body terms; it was success- 
fully employed to study the properties of silicon [1, Q 
(just to name a few), and has recently been used as a 
coarse-grained model of water In this letter, we 

examine this system and identify the different configu- 
rations that can occur around a particle (its first coor- 
dination shell). We then assign energies, volumes and 
degeneracies to each possible configuration, or "species" , 
using simple theoretical considerations. Finally, we cal- 
culate the free-energy of an ensemble of these species un- 
der constant volume and temperature. Since our calcula- 



tion assumes constant volume and not constant pressure, 
it is natural to consider experiments were the liquid is 
confined to nano-pores ll|. We demonstrate that such 
a theory, based on a statistical mechanical description of 
the local structure of the liquid can explain the anoma- 
lous behavior of tetrahedral liquids due to a competi- 
tion between energy and entropy. We also explain why 
the same mechanism that is responsible for the anoma- 
lous behavior can, in principle, create a phase-transition 
when certain conditions are met. An important feature of 
tetrahedral liquids, as observed in simulations and exper- 
iments, is that the average coordination numbers changes 
significantly as a function of the temperature. For that 
reason we would like to propose an approximate statis- 
tical mechanics for the Stillinger- Weber liquid, based on 
a description of the local structures. This kind of theory 
was proposed a long time ago for liquids by Bernal and 
others 12 , 13| , and recently it was developed to describe 
the structural changes that occur in supercooled liquids 
14] . This kind of "upscaling" is based on an assumption 
of the additivity of the energy of the local structure and 
on the fact that supercooled liquids are ergodic |15| . 

The Stillinger- Weber model Q is a widely used em- 
pirical model of silicon which assumes that the leading 
terms in the N-body quantum mechanical interaction are 
the two-body and three-body terms. The energy function 
is: 
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where is a Lennard- Jones-like potential with a finite 
cutoff, and v-$ is constructed such that it nullifies when 
the angles between the three neighbors are 9 = 109.47°. 
This property guarantees that the system has a diamond- 
cubic ground-state at atmospheric pressure. The exact 
details of the functions V2 and v% can be found in ref. Q . 
In order to take into account the contribution of the three 
body interactions we make the following approximation: 
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we assume that the first shell of neighbors of each par- 
ticle can be mapped to the first-shell of a particle in a 
BCC lattice oriented in an arbitrary direction. This is an 
approximation since in the liquid, particles can assume 
arbitrary positions. The BCC first-shell seems to be the 
best discretization of the liquid first-shell environment 
since the central particle has at most 8 neighbors. This 
is consistent with simulations that use the SW potential, 
since it is possible to arrange 4 particles around the cen- 
tral particle in a way that replicates the environment of 
a particle in a diamond-cubic lattice (the ground-state of 
the Stillinger- Weber potential) . The values of the ener- 
gies of the different states i were calculated from the SW 
potential by placing particles at the different possible lo- 
cations in the BCC first-shell. The energies were found 
to be: 

Ei= {0, -e , -2e , -2e + ei, -3e , -3e + 2ei, 
-4e , -Ae Q + 3d, -4e + 4e x , -5e + 4e x , 
— 5eo + 6ei, -6e + 8ei, — 6eo + 9ei, 
-7e + 12ei, -8e + 16ei}, (2) 

where eo is the typical energy scale of the two-body in- 




FIG. 1: Color online: representative nearest neighbor config- 
uration for n = 2 (a) and three possible configurations for 
n = 4 (b)-(d), each having different three-body contributions 
to the potential energy. The central particle is in green, the 
neighbors are colored red and the rest of the possible particle 
positions are semi-transparent. 

teraction, and ei is the energy scale of the three-body 
term. At this point eo and ei are left as free parameters. 
There are several different configurations which give the 
same energy and thus are degenerate. For a particle with 
no nearest neighbors there is only one possible configura- 
tion and thus the degeneracy of the state is go = 1. When 
we have one nearest-neighbor in the first-shell, there are 
eight ways to arrange the neighbor in the BCC structure 
and gi — 8. When there are two nearest-neighbors there 
are twelve ways to arrange the particles in a configura- 
tion where the three-body term does not contribute to 



the energy, and sixteen different ways to arrange the two 
neighbors in a way were there is one contributions of the 
three-body term (the three-body term is non-zero only 
when there are at least three particles involved in the in- 
teraction). Similarly, we calculated the degeneracies that 
corresponds to the energies in Eq. ^ : 

gi = {1, 8, 12, 16, 8, 48, 2, 32, 36, 8, 48, 12, 16, 8, 1} . (3) 

Each species occupies a typical volume in space, which 
depends on the geometry of the local configuration and 
the forces which the neighboring particles apply on the 
central particle. The typical volume per species was cal- 
culated by minimizing the energy for the different config- 
urations and calculating the volume of the Voronoi cell 
of each particle using the voro++ package given in [lij ]. 
The volume of each species increases with temperature 
(with a rate a) due to thermal fluctuations while the 
three-body interactions cause the species to decrease its 
volume with temperature at a rate (3 that depends on 
the number of three-body terms in the interaction. This 
is due to the repulsive nature of the three-body term 
which is more important at lower temperatures (fluctua- 
tions around the angle of minimum energy contribution 
are smaller). Both of these contributions were combined 
to give: 

Vi (T) = { 1.395 +aT, 1.297+ aT, 1.199 + aT, (4) 
1.220 — f3T + aT, 

1.101 +aT,l.l4Q-20T + aT,l.O + aT, 
1.058 -3/3T + aT, 1.083 -4/3T + aT, 
0.974 -3/3T+aT, 1.019- 6/3T+ aT, 
0.943 -8/3T + aT,0.968- 9/3T + cvT, 
0.910 -12 + 0.870- ISfiT + aT}, 

where a and f3 have units of [if -1 ]. The volumes are nor- 
malized by the volume of a particle in the ground-state at 
zero pressure (four neighbors in the Diamond-Cubic lat- 
tice - no contribution due to three-body interactions) and 
are dimensionless. The parameters a and /? together with 
eo and t\ complete the parametrization of the model. The 
statistical ensemble of particles at constant volume and 
temperature is given by the Helmholtz free energy: 

14 14 N 

F = U -TS = y^N k E k +Ty^N k ln—±-, (5) 

where S — — J2k~o^ k ^ n (^ k /^9 k ) ^ s tne configura- 
tional entropy. The free energy will have a minimum for 
the equilibrium values of the concentrations of species 
Pk = -^r- In order to account for the isochoric condition 
(constant volume) we need to demand that the total vol- 
ume equals V. We do that by adding a volume constraint 
together with a Lagrange multiplier A: 

14 14 14 

/= ^Pk E k +T^2p k \n— + \(^PkVk -v), (6) 

k=0 k=0 ® k fc=0 
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where v = is the average volume per particle and / = 
F/N is the free-energy per particle. After minimizing 
the free-energy under the constraint we get the modified 
Boltzmann distribution: 

Pi(T) = ^ 9l e-^ , Z{T) = ^> e -T, (7) 

i 

where Z{T) is the partition function. To that we need 
to add the constraint: 



(8) 



Note that different states i may have the same coordi- 
nation number j. We define the concentration Cj(T) of 
the species having coordination number j at temperature 
T . Solving numerically Eq. ([7]) under the constraint ((5J) 
we can find the average concentrations of the different 
species as a function of the temperature, from which we 
can compute c, (T) . To asses the quality of the theory we 
present results from the NVT molecular dynamics simu- 
lations which were averaged over fifty different ensembles 
each containing 4096 particles, see Fig. [5J The theoreti- 
cal predictions of the model, are obtained by fitting the 
four free parameters (eo, e±, a and 0). The results of 
are shown in the inset of Fig. [2] Obviously the the- 
ory provides a semi-quantitative description of the nu- 
merical data. Once we have fitted the free parameters 
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FIG. 2: Color online: Concentrations Cj(T) of particles with 
coordination numbers j from simulations. Inset: the theoret- 
ical results for comparison. The fitted values of the parame- 
ters are e = 8.7 • 10" 6 , ei = 3.1 ■ 10" 6 , a = 1.4 • 10" 5 and 
=4.2- 10" 6 . 

to agree with the plots of the concentrations as a func- 
tion of temperature, we do not change them, and look 
at the model as fixed to assess its thermodynamic prop- 
erties. As stated above, tetrahedral liquids exhibit pe- 
culiar low temperature behavior, such as a density max- 
imum, an increase in the specific heat capacity and a 



negative thermal expansion. In figure [3^, we observe the 
calculated density p for P = — ^ =0. A maximum 
appears at about T = 950-fT, roughly the same temper- 



ature as observed in 17J for molecular dynamics of the 
SW model. We also observe a small density minimum 
at around T = 700-ftT which is in agreement with 11711 
and with a recently published experiment on water |18| . 
The thermal expansion and specific heat capacity at zero 
pressure are also shown in figures^ and[3fc. The specific 
heat shows a large increase while the thermal expansion 
decreases and becomes negative, as observed experimen- 
tally in water. 
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FIG. 3: Density maximum and minimum, thermal expansion 
minimum and specific heat maximum for P = 0. 

To understand these results, we consider the energy 
and the entropy as a function of the volume V at two dif- 
ferent temperatures (T = 755K and T = 30187^). These 
can be immediately computed from the model, and pro- 
vide us with insights into the physical mechanism that 
allows for the anomalous behavior shown in the previ- 
ous subsection. Observe Fig. \Q: at high temperatures, 
the configurational entropy (divided by 10 for compari- 
son) has one maximum at a relatively high volume. The 
reason is that at a higher volume there are more ways 
to arrange species with different volumes. Since at high 
temperatures the term — T S in the free energy is dom- 
inant, the minimum of the free-energy will appear at a 
higher volume. Lowering the temperature means that 
the potential energy becomes more important than the 
entropy and the equilibrium configurations tends to pre- 
fer the grounds-state species which have the lowest en- 
ergy. Therefore, when the temperature is low the energy 
will prefer the volume V=l and will develop a minimum 
around that point, as one can observe at Fig. [4j\ When 
we constrain the ensemble to have a higher volume, more 
species appear in order to fill the gaps in space; simi- 
larly smaller species appear when the volume is smaller 
than unity (cf. Fig. The configurational entropy also 
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develops a minimum around V = 1 since most of the 
particles at this volume will only have four neighbors. 
Therefore, the minimum in the free-energy decreases and 
approaches the ground-state volume as the temperature 
is lowered. As we explained above, the increase in the 
concentrations of the various species when the volume is 
altered with respect to V = 1, increases the configura- 
tional entropy which causes it to exhibit a double max- 
imum, where a small increase in the volume results in 
a large increase in the entropy (cf. Fig. 0J-). If this 
effect were sufficiently strong, the entropy could "pull" 
the free-energy to higher volumes, despite the low tem- 
perature. This causes the density to have a maximum 
at a low temperature. If the increase in entropy due to 
volume changes were even larger, the free energy could 
develop a second minimum - a first order phase transi- 
tion. In the present model the entropy contribution is 
not large enough to cause a phase transition. Beaucage 



at al. 17] found that the liquid- liquid phase-transition 



in the Stillinger- Weber model is very sensitive to small 
perturbations of the potential and simulation parame- 
ters. When the three-body contribution increases by as 
little as 5%, the phase transition disappears while the 
system maintains the density-maximum property. Their 
result indicates that similar potentials can lead to a phase 
transition in some cases but not in others. A liquid can 
exhibit anomalous behavior without undergoing a phase 
transition. We suggest that this insight applies to real 
tctrahcdral liquids: some may exhibit a phase transition 
while others may show anomalous behavior without a 
phase transition. 
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FIG. 4: Color online: (a) Concentrations of Cj(T) as a func- 
tion of volume for T = 755K. (b) Comparison of the potential 
energy and the entropy for two different temperatures. 

In summary, we have constructed a simple statistical 
mechanics theory which describes the local structure of a 
generic model of tetrahedral liquids, the Stillinger- Weber 
model. We showed that the theory, intended originally to 
describe only the structure, can account for the anoma- 
lous behavior of tetrahedral liquids and that this behavior 
can be understood as a competition between configura- 



tional entropy and potential energy. We suggest that the 
same physics can give rise to a liquid-liquid phase transi- 
tion, and that it is possible that in some tetrahedral liq- 
uids the anomalous behavior is accompanied by a phase- 
transition while in others the phase transition does not 
occur. It is important to stress however that the existence 
of a phase transition is not a pre-requisite for the ther- 
modynamic anomalies which are perfectly present even 
without any phase transition in the supercooled liquid. 
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